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Broken Symmetries, Zero-Energy Modes, and Quantum Transport in Disordered 
Graphene: From Supermetallic to Insulating Regimes 
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The role of defect- induced zero-energy modes on charge transport in graphene is investigated using 
Kubo and Landauer transport calculations. By tuning the density of random distributions of mono- 
vacancies either equally populating the two sublattices or exclusively located on a single sublattice, 
all conduction regimes are covered from direct tunneling through evanescent modes to mesoscopic 
transport in bulk disordered graphene. Depending on the transport measurement geometry, de- 
fect density, and broken sublattice-symmetry, the Dirac point conductivity is either exceptionally 
robust against disorder {supermetallic state) or suppressed through a gap opening or by algebraic 
localization of zero-energy modes, whereas weak localization and the Anderson insulating regime 
are obtained for higher energies. These findings clarify the contribution of zero-energy modes to 
transport at the Dirac point, hitherto controversial. 



The electronic transport properties of graphene are 
known to be very peculiar with unprecedented manifes- 
tations of quantum phenomena such as Klein tunneling 
[iIjI^I, weak antilocalization [3|,|4|, or the anomalous quan- 
tum Hall effect [5|, |6| , all driven by a 7r-Berry phase stem- 
ming from graphene sublattice symmetry and pseudospin 
degree of freedom [3-Q- These fascinating properties, 
yielding high charge mobility [10, \lM j ^.re robust as long 
as disorder preserves a long range character. The funda- 
mental nature of transport precisely at the Dirac point 
is, however, currently a subject of fierce debate and con- 
troversies. Indeed, for graphene deposited on oxide sub- 
strates, the nature of low-energy transport physics (as its 
sensitivity to weak disorder) is masked by the formation 
of electron- hole puddles [9|. A remarkable experiment 
has, however, recently demonstrated the possibility to 
screen out these detrimental effects 1^, providing ac- 
cess to the zero-energy Dirac physics. An unexpectedly 
large increase of the resistivity at the Dirac point was 
tentatively related to the Anderson localization [l2|, |l3[ 
of an unknown physical origin and questioned interpre- 
tation [3]. 

Of paramount importance are therefore the low-energy 
impurity states known as zero-energy modes (ZEMs) 
[iSl |l6[, whose impact on the Dirac-point transport 
physics needs to be clarified. ZEMs are predicted or 
observed for a variety of disorder classes, as topologi- 
cal defects (mainly vacancies) lla. llTII. adatoms cova- 
lently bonded to carbon atoms [18, llSf, and extended 
defects as grain boundaries [20, l21|. As recently con- 
firmed by scanning tunneling microscopy experiments on 
graphene monovacancies |22| . ZEMs manifest as wave 
functions that decay as the inverse of the distance from 
the vacancy, exhibiting a puzzling quasilocalized charac- 
ter, whose consequences on quantum transport remain. 



to date, highly controversial. First, ZEMs have been pre- 
dicted to produce a supermetallic regime by enhancing 
the Dirac-point conductivity above its minimum ballistic 
value CTinin = ^e^ /irh [23|,|2J|, an unprecedented conduct- 
ing state, which could be, in principle, explored experi- 
mentally |25l427| . Second, a similar increase of the Dirac- 
point conductivity with the defect density has been also 
reported in the diffusive regime of two-dimensional disor- 
deredgraphene in the presence of vacancies or adatoms 
|l9l l28j . These results contrast with the semiclassical 
conductivity found with the Boltzmann approach ,171,291- 
|32|, and suggest the absence of quantum interferences 
and localization effects observed for other types of disor- 
der 33-35]. Finally, transport experiments in intention- 
ally damaged graphene also report on puzzling conductiv- 
ity fingerprints, whose physical origin remains to be fully 
understood 36|, |37| . A comprehensive picture of the role 
of ZEMs on quantum transport properties in disordered 
graphene is therefore crucially missing and demands for 
further theoretical and experimental inspection. 

This Letter provides an extensive analysis of the con- 
tribution of zero-energy modes to quantum conduction 
close to the Dirac point in disordered graphene. Us- 
ing the Kubo-Greenwood and Landauer transport ap- 
proaches, different regimes are numerically explored by 
changing the aspect ratio of the transport measurement 
geometry, and by tuning vacancy density and sublattice 
symmetry breaking features. The robustness of the su- 
permetallic state induced by ZEMs is shown to be re- 
stricted to very low densities of compensated vacancies 
(equally distributed among both sublattices). This oc- 
curs as long as tunneling through evanescent modes pre- 
vails. In the absence of contact effects, an increase of 
the conductivity above Ae'^/nh is obtained for the semi- 
classical conductivity at the Dirac point and ascribed to 



a high density of ZEMs, but the quantum conductivity 
analysis unequivocahy reveals a localization regime. For 
a totally uncompensated vacancy distribution (populat- 
ing a single sublattice), the delocalization of ZEMs in 
real space is strongly prohibited for a large energy win- 
dow around the Dirac point owing to the formation of 
a gap, whereas no appreciable difference of high energy 
transport (above the gap) is found compared with the 
compensated vacancy case. We would like to mention 
that some interesting cases of uncompensated impurities 
and defects have been reported experimentally [38|-|40J, 
whose results demand further exploration. 

System description and methodology.- We consider a 
finite concentration n of vacancies either distributed at 
random exclusively on one of the two sublattices {ua = n, 
the number of vacancies per carbon atoms in sublattice 
A and ub — 0, uncompensated case), or equally dis- 
tributed vacancies on both sublattices (ua — ub — fi/2, 
compensated case). The electronic and transport prop- 
erties are investigated by using a tight-binding model 
with a single p^ orbital per atom and first nearest neigh- 
bor coupling. We model the vacancies by removing the 
corresponding orbitals from the Hamiltonian 
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To investigate the various transport regimes, two com- 
plementary approaches are used. For studying two- 
dimensional (bulk) disordered graphene, real-space quan- 
tum wave packet dynamics and Kubo conductivity are 
calculated [33|, |3J, l4ll444j . The zero- frequency conductiv- 
ity cr(-E, t) for energy E and time t is given by a{E, t) ~ 
e'^p{E)AX'^{E,t)/t, where p{E) is the density of states 
(DOS) and AX'^{E,t) is the mean quadratic displace- 
ment of the wave packet at energy E and time t: 



AX^{E,t) 



TT[5{E-n)\Xit)-X{0)f' 
TiiS{E - H)] 



(1) 



A key quantity in the analysis of the transport properties 
is the diffusion coefhcient D{E, t) = AX'^{E, t)/t. In dis- 
ordered systems, D{t) generally starts with a short-time 
ballistic motion followed by a saturation regime, which al- 
lows us to estimate the transport (elastic) mean-free path 
£e from its maximum value as (.e{E) = D'^'^^{E)/2v{E) 
where v{E) is the velocity. The semiclassical conductiv- 
ity cTsc is given analogously by the maximum conductiv- 
ity. Depending on disorder strength, D{E,t) is found to 
decay at longer times owing to quantum interferences, 
whose strength may dictate weak or strong (Anderson) 
localization at the considered time scale. Calculations are 
performed for systems containing several millions of car- 
bon atoms, allowing the capture of all relevant transport 
regimes. We also study the ballistic limit of transport 
through finite graphene samples, by considering strip ge- 
ometries with width W and length L (with W/L ^ 1) 
between two highly doped semi- infinite ribbons (of identi- 
cal width). This two-terminal transport geometry gives 
access to the contribution of ZEMs in graphene trans- 
port when the charge flow is conveyed by contact-induced 
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FIG. 1. Main frame: 
(compensated case): 



Conductivity of graphene with n—0.8% 
semiclassical value asc (solid line). 



o"min = 4e /nh (dotted line), and Kubo conductivity at vari- 
ous time scales. Left inset: DOS for varying vacancy density, 
together with the pristine case (dashed line). Right inset: 
Mean-free paths for n = 0.1%, 0.2%, and 0.4%. 



evanescent modes. The doping of contacts is simulated 
by adding an on site energy of -1.5 eV to the corre- 
sponding orbitals, which generates a large DOS imbal- 
ance between the contacts and the central strip at the 
Dirac point (E = 0). The zero-temperature conductiv- 
ity of the graphene strip is then computed as <7{E) = 
{2e'^/h)T{E)L/W, where T{E) is the transmission co- 
efficient evaluated within the Green's function approach 
[45l |46| . When L <C W, low-energy transport is domi- 
nated by tunneling through the undoped region yielding 
a universal ballistic value a{E ~_0) ~ CTm n — ie'^/irh at 



the Dirac point for clean strips |25l - l27l |45 1 . 

ZEMs effects in two-dimensional disordered graphene.- 
We start by considering the compensated case, which 
globally preserves the sublattice symmetry. Figure [T] 
(left inset) gives the density of states of the system as 
a function of the energy E for different vacancy densi- 
ties 71. In agreement with prior results [15|, llfll, the DOS 
shows the rise of a broad peak around E = 0, which 
witnesses the presence of ZEMs generated by disorder. 
Their nature, however, is not encoded in this feature but 
needs to be analyzed by studying transport characteris- 
tics such as the mean-free path (Figll] right inset) and 
conductivity (Figll] main frame). The mean- free path (.^ 
is seen to be strongly energy dependent with minimum 
values close to the Dirac point, as expected for short- 
range scatterers [43|, |4j| . By increasing the vacancy den- 
sity within the range [0.1%, 0.4%], (.^ drops from tens of 
nanometers down to few nanometers, and roughly varies 
as ^e ~ l/?^! which is in agreement with the Fermi golden 
rule. Interestingly, we find for the semiclassical conduc- 
tivity Gsc ^ E for a high enough energy (above 0.3 eV 
for n = 0.8%), whereas it saturates to (Tmin at low en- 
ergy with a higher value around the Dirac point owing to 




FIG. 2. Main frame: (Tsc{E) and a{E,t) for graphene (uncom- 
pensated case) and energy resolution 77 = 3 meV. Left inset: 
DOS with an energy gap revealed by rj scaling [52!] and ZEMs. 
Right inset: Diffusion coefficients at E — 0.5 eV and E = 
{rj = 3 meV) for both compensated (AB) and uncompensated 
(AA) cases. All data for n = 0.8% 



the DOS enhancement induced by midgap states. When 
increasing the vacancy density, the minimum conductiv- 
ity Ae^ /irh around the Dirac point extends over a larger 
energy region (not shown here). 

The obtained short £e and minimum semiclassical con- 
ductivities suggest a strong contribution of quantum in- 
terferences, which is further evidenced by the decay of 
the Kubo conductivity below (Jmin for sufficiently long 
time scales, see FigH] (main frame). Depending on the 
energy, the observed downscaling of the quantum con- 
ductivity versus time can be described by a logarith- 
mic correction (weak localization), an exponential de- 
cay (strong localization), or by algebraic localization of 
the ZEMs. As detailed in the Supplemental Material, 
the quantum correction to the conductivity [Sa{X) ~ 
a{X)—(Jsc] at E=OA eV is numerically found to downscale 
as 5a{X) - -2e^/{7rh) ln(A/Ae) [with A ee ^/AX^{t) the 
time-dependent wave packet space extension and Ae re- 
lated to (e [13] • Differently, for £^=0.2 eV the length- 
dependent conductivity exhibits an exponential behavior 
a ^ exp(— A/f) (^ the localization length), evidencing 
a strong- localization regime 13]. Exactly at the Dirac 
point, the conductivity decays following a power law 
a ^ X~^ . This behavior is actually in full agreement 
with the localization of the ZEMs measured experimen- 
tally by scanning tunneling microscopy 22]. 

A remarkably different picture emerges in the uncom- 
pensated case, for which the sublattice symmetry is fully 
broken. The DOS shown in Figl2] (left inset) evidences 
the presence of ZEMs sharply peaked at E = 0. fn 
contrast to the compensated case, the depletion of the 
low-energy conductivity is here inherited from the pres- 
ence of energy gaps lajllSl. The semiclassical conductiv- 
ity strongly increases when approaching the Dirac point. 
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FIG. 3. Main frame: Conductivity for strips with W — 150 
nm, L — 15 nm, and a compensated vacancy density up to 
2%. Inset: Same information for uncompensated vacancies 
with densities up to 1%. 



much more than in the compensated case and also in- 
creases when improving the energy resolution. However, 
the large value of age does not reflect the extendedness 
of the corresponding ZEMs. This can be rationalized by 
scrutinizing a{E = 0,t) and D{E = 0,i), which are ac- 
tually strongly decaying with time. Indeed D{E = 0,i) 
becomes extremely small compared to that at finite en- 
ergies (e.g. at 0.5 eV) and much smaller compared to 
the compensated case with the same vacancy concentra- 
tion (see FiglH right inset). Additionally, D{E = 0,i) 
decays when improving the energy resolution (not shown 
here) , thus demonstrating that although many ZEMs are 
present, they do not participate in conduction, and that 
the large value of cr^c obtained numerically results from 
the high DOS at E ~ 0. Furthermore, the physical rel- 
evance of a semiclassical conductivity at the Dirac point 
is highly questionable (see also Supplemental Material). 
For the quantum conductivity, on the other hand, the 
strong decay of a{E — 0,t) with time is consistent with 
localized modes similar to the compensated case. We also 
find that away from the Dirac point a higher energy reso- 
lution reduces ctsc and a{t) as observed for the DOS, thus 
unambiguously indicating the energy gap as the origin of 
the conductivity decrease, and ruling out any diffusive 
regime and Anderson localization phenomenon. Finally, 
for larger energies away from the gap region, one observes 
that the wave packet dynamics for the compensated (AB) 
and uncompensated (AA) cases are very similar; see Figl2] 
(right inset). This discards any singular transport mech- 
anism in an uncompensated situation, which is different 
from previous reports on hydrogenated graphene 



ZEMs effects in disordered finite graphene strips.- In 
contrast to two-dimensional graphene, the role played by 
ZEMs in transport through finite strips in between highly 
doped contacts turns out to be quite different. In this 
configuration, the contacts have much a higher density 
of propagating states than the central strip, especially at 
the Dirac point. Accordingly, many states from contacts 



tunnel through the strip as evanescent modes, yielding 
a minimum ballistic value CTmin = 4e^/7r/i for clean sam- 
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pies J25l - l27l |. The presence of ZEMs increases the num- 
ber of available states at the Dirac point in the central 
strip. Two competing transport mechanisms then drive 
the conductivity behavior, namely an enhanced tunnel- 
ing probability assisted by ZEMs together with multiple 
scattering and quantum interferences, which develop ow- 
ing to the randomness of vacancies distribution. 

Figure[3](main frame) shows the quantum conductivity 
a for a strip with length L = 15 nm, width W = 150 nm, 
and compensated vacancy density in the range [0%, 2%]. 
In the absence of vacancies, a shows the minimum con- 
ductivity a{E = 0) = (To ~ Cmin expected for the ballistic 
limit when L <^ W (see the horizontal dotted line) [25j. 
For n = 0.1%, the strip length is close to the mean- free 
path; see FiglTJ Therefore, the transport along the strip 
remains quasiballistic, a fact further confirmed by the 
smooth decay of a all over the spectrum except at the 
Dirac point, where cr keeps a larger value. For higher 
densities and away from the Dirac point, the decay of 
a{E) with n is consistent first with the occurrence of 
a diffusive regime and then with localization phenom- 
ena, as revealed by the strongly fiuctuating conductivity. 
Note that despite the few nanometers short mean-free 
path, even for n = 2%, the conductivity remains signifi- 
cant as a consequence of the large number of conductive 
channels that penetrate the undoped strip. The con- 
ductivity around the Dirac point is further scrutinized in 
FigH] (bottom inset) for strips with L = 15 nm, W = 150 
nm, and compensated vacancy densities up to n = 1%. 
To reduce sample-to-sample fiuctuations, all the results 
were averaged over 20 random disordered configurations. 
Far from the Dirac point, the conductivity is found to 
decrease regularly with n. At E = 0, notably enough, 
a peak is always present, which can slightly exceed ctq 
at very low density {n ^ 0.04%). This indicates that 
the ZEMs generated at the Dirac point are suSiciently 
delocalized to assist (and even enhance) electron tunnel- 
ing through the strip. Backscattering becomes eventu- 
ally dominant for sufficiently high defect concentration, 
as manifested by the smooth conductivity decrease. The 
dependence of the conductivity peak (cpoak) on the differ- 
ent system parameters is reported in Fig|4] (main frame) 
for compensated vacancy densities up to 5% and lengths 
up to 15 nm. The decrease of (Tpeak with n is very slow, 
especially for the shortest strip, and even for strong dis- 
order (n = 5%) CTpcak remains significantly large. As il- 
lustrated in Fig|4](top inset), Cpeak is actually a universal 
function oi nx L^. Remarkably enough, CTpoak fluctuates 
around or goes slightly above ao for very low nxL? < 10, 
thus supporting the possibility for a supermetallic state, 
introduced by Ostrovsky and co-workers [23|, |2J]. For 
n X L^ ^ 10, (Tpoak decreases roughly logarithmically, as 
the result of finite size effects and proximity between va- 
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FIG. 4. Main frame: Average conductivity peak versus n 
for strips with W = 150 nm and L — 5, 10, and 15 nm. 
The shaded areas around the curves indicate the standard 
deviation with respect to the average value. Top inset: Same 
as the main frame but as a function of n x L^. The thick 
straight hue is a guide to the eye. Bottom inset: Average 
conductivity for W — 150 nm, L = 15 nm, and various n. 



The conductivity of graphene strips (with M^ = 150 
nm, L =15 nm, and n up to 1%) for uncompensated 
vacancies is reported in Fig|3] (inset). In marked con- 
trast with the prior case, a gap develops at low density 
together with a reduced but finite conductivity peak at 
E — 0. As for the case of two-dimensional graphene 
(Figl2), the gap formation leads to the suppression of 
tunneling due to the almost vanishing DOS. The Dirac 
conductivity peak is a signature of the highly localized 
nature of zero-energy states generated by uncompensated 
vacancies [161, which are not enough spatially extended to 
significantly contribute to tunneling and obviate to the 
DOS decrease. More details on the energy gap scaling 
and, in general, on the transport properties away from 
the Dirac point will be published elsewhere [49J. 

In conclusion, the contribution of ZEMs to quan- 
tum transport in disordered graphene has been dis- 
cussed for various transport geometries and sublattice 
symmetry-breaking situations. Our findings provide a 
broad overview of the low-energy transport phenomena 
in graphene in the presence of ZEMs, including the for- 
mation of an insulating state at the Dirac point, accessi- 
ble in absence of electron- hole puddles [l^j- The role of 
electron-electron interaction (here neglected) , might also 
play some important role in capturing the full picture 
and deserves further investigation 50|, |5l| . 
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SUPPLEMENTAL MATERIAL 

BROKEN SYMMETRIES, ZERO-ENERGY MODES AND QUANTUM TRANSPORT IN DISORDERED 

GRAPHENE: FROM SUPERMETALLIC TO INSULATING REGIMES 

by Alessandro Cresti, Frank Ortmann, Thibaud Louvet, Dinh Van Tuan, and Stephan Roche 



In this Supplemental Material, we demonstrate that 
the decay of the quantum conductivity shown in Fig. 
1 (main frame) of the main paper for compensated va- 
cancy concentration of 0.8% (at finite energies) follows 
the conventional scaling theory of localization for two- 
dimensional disordered systems. To this end, Fig. 31] 
shows the conductivity as a function of the wave-packet 
spread A(i) — y^AX^{t) at different energies and for dif- 
ferent energy resolution parameters rj. 

In Fig. HI](a), we see that the decay of the conductivity 
with A observed at £^ = 0.4 eV (due to the quantum 
correction) can be fitted by — l/7rln(A/Ae) (using units 
2e^/h and with Ae related to the mean free path [1[). 
This clearly indicates that the system is in the weak- 
localization transport regime. In the same panel (a) a 
strongly reduced broadening 77 = 0.8 meV yields the same 
conductivity. 

At lower energy {E — 0.2 eV), a strong localization 
regime is obtained as seen in Fig. 31] (b), for different 
energy resolutions from rj — 5 meV down to rj = 0.4 
meV. The length-dependent conductivity decays expo- 
nentially as cr ~ exp(— A/^), hence evidencing the strong- 
localization regime (^ the localization length). Note that 
this scaling law is observed independently of the energy 
precision parameter 77, thus indicating that our approach 
is able to unambiguously catch the physics of the system 
and that there is only a residual quantitative, but not 
qualitative, dependence on 77. Moreover, the localization 
length varies only weakly at lowest 77 indicating a limit 
^ wlO nm when 77 — > 0, which confirms the reliability of 
the numerical simulation. 

At the Dirac point {E = 0), localization is observed 
in Fig. 31] (c) since the conductivity decays with length 
A. However, in contrast to finite energies, it follows a 
power-law behavior a oc A" with a < 0. The inset shows 
a upon decreasing the broadening 77 down to about 0.4 
meV, which is the present limit of our numerical resolu- 
tion. Note that the observed behavior is consistent with 
the limit a = —2 for 77 — >■ 0, which has been observed 
experimentally [3] for the localization of ZEMs by means 
of scanning tunneling spectroscopy. The localization at 
-E = is even stronger than in the Anderson regime and 
can therefore not be attributed to multiple scattering and 
quantum interference effects, i.e. the strong localization 
regime, but is rather a signature of zero-energy modes. 
This is further corroborated by the length A '--^ 5 nm 
over which a localizes, which is on the same order as 
the spatial extension of the bound states experimentally 
measured ^. 

We would like to point out here that our results for 
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FIG. SI. Length-dependent conductivity for different ener- 
gies and 0.8% vacancy concentration in the compensated case. 
(a) Conductivity a and quantum correction Sa = a — asc at 
E — 0.4 eV. The logarithmic fit confirms the weak-localization 
regime, (b) Low energy conductivity {E = 0.2 eV) and cor- 
responding fit indicate Anderson locafization regime, (c) At 
zero energy the conductivity decay is even stronger and can- 
not be fitted with an exponential decay, (d) Conductivity at 
largest simulated times (8.2ps) and its residual dependence 
on 7;. 



compensated vacancies are well-defined and converge in 
the limit of small rj. Figure 31] (d) finally shows that at 
the largest time considered for the calculation of the con- 
ductivity (8.2ps), (j{E) is well controlled when decreasing 
77, with a more pronounced noise level at smaller 77, an ef- 
fect which defines a lower limit for 77 to avoid non-physical 
mathematical singularities. 

The case of uncompensated vacancies is more subtle 
and the role of 77 depends on the region of the energy 



spectrum considered. At high energies, as seen in Fig. 2 
of the main paper, the physics is the same as for com- 
pensated vacancies and the hmit 77 — > leads to quan- 
titatively robust results in the weak-localization regime. 
At lower energies, but not at the Dirac point, the sys- 
tem exhibits a genuine energy gap where electrons cannot 
propagate, thus leading to vanishing conductivity, which 
is therefore not related to a transition from diffusive to 
localized regime. The conductivity suppression is there- 
fore of different nature as compared to the compensated 



to 0.4 meV) and with a steeper decay for smaller 77, which 
is consistent with bound states of vanishing conductivity 
similar to the compensated case. In contrast, we observe 
that the semiclassical conductivity diverges with small 
rj. The reason is that, for the uncompensated case, all 
vacancy-induced modes are exactly at -E = and their 
corresponding DOS and semiclassical conductivity have a 
(5-like distribution centered in the gap where no propaga- 
tion is possible. However, the broadening and the height 
of the DOS peak (as well as Usc peak) are artificially 
driven by the finite parameter ry. 



At the Dirac point, our numerical approach provides an 
ry-dependent value of the quantum conductivity that pro- 
gressively decreases with 77, consistently with vanishing 
conductivity. We find that, similarly to Fig. 311(c), the 
conductivity decays on a very short length scale, which is 
even shorter than for the compensated case. This is ob- 
served for all the ry values used in the simulations (down 
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